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SUMMARY 


A numerical method will be presented to solve the three-dimensional Navier-Stokes equations 
in combination with a full Reynolds-stress turbulence model. Computations will be shown for three 
complex flowfields. The results of the Reynolds-stress model will be compared with those predicted 
by two different versions of the k - u model. It will be shown that an improved version of the k - uj 
model gives as accurate results as the Reynolds-stress model. 


Faster computer performance and improved numerical methods not only allow the numerical com- 
putation of more complex flowfields, but also enable the use of higher order turbulence models in 
Reynolds averaged Navier-Stokes computations. This is reflected by an increasing number of industrial 
applications utilizing one-equation or two-equation models instead of the much cheaper algebraic models. 
Ihese improvements allow the research community to go one step further and to adopt full second- 
oidei closuie models (full Reynolds-stress models) for the computation of complex three-dimensional 
flowfields. 

The first part of the paper deals with the implementation of a Reynolds-stress model (RSM) into 
a Navier-Stokes code. The Navier-Stokes method is the INS 3D code developed by S. E. Rogers and 
D. Kwak (ref. 1). It solves the Reynolds-averaged incompressible Navier-Stokes equations in general 
curvilinear coordinates with a pseudo-compressibility method. The RSM used in this study is the 
multiscale (MS) model of D. C. Wilcox (ref. 2). In this model the Reynolds-stresses are partitioned 
into two sets of scales. The large eddies, which act inviscidly, and the small eddies that are isotropic 
and provide the dissipation. The main advantage of this partitioning is that no diffusion terms appear 
in the equations for the large eddies, at the expense that one more equation has to be solved, compared 
to a single scale RSM. It should be noted that the most important part of the model, the pressure-strain 
term, is identical to the Launder-Reece and Rodi model (ref. 3). 

The implementation of an RSM into a Navier-Stokes code is not a trivial matter and can lead to 
severe stability problems (refs. 4 and 5). The discretization and the linearization of the equations is 
of crucial importance to the success of the computations. In the present approach, the equations are 
discietized with a higher order upwind scheme. In order to enhance stability, a total variation diminishing 
(TVD) (ref. 6) limiter is imposed on the numerical fluxes. Furthermore, the eddy-viscosity formulation 
of the underlying two-equation model is retained in the implicit side of the mean flow equations to 
enhance diagonal dominance. This approach is similar to that of an apparent eddy-viscosity adopted in 
lefeience 4, except that the apparent viscosity is isentropic in the present discretization. 

In the second part of the paper, results will be shown for four increasingly complex flowfields. The 
first flow is a zero pressure gradient flat plate boundary layer, which is mainly included as a consistency 
check for the numerical implementation. The second problem is an adverse pressure gradient flow with 
separation (refs. 7 and 8). The third is a flow over a circular cylinder, with a rotating and a stationary 
section (refs. 8 and 9) and the last is a fully three-dimensional pressure-driven boundary layer flow 
(ref. 10). All flowfields have also been computed with the k-uj two equation model in its original form 
and in a modified veision recently proposed by the author for adverse pressure gradient flows (ref. 11). 



Each of the flowfields focuses on a different aspect of the potential superiority of the RSM over 
the two-equation models. In the adverse pressure gradient flow, the main difference between the k — to 
model and the RSM is that the latter accounts properly for' the transport of the Reynolds stresses. The 
spinning cylinder flow will investigate the potential of the RSM to predict strong curvature effects 
that are known to be beyond the scope of eddy-viscosity models. The three-dimensional pressure- 
driven boundary layer flow will show whether the RSM is able to properly predict the differences 
between the stress and the strain direction as well as the reduction of the turbulent shear-stress observed 
experimentally in three-dimensional flows (refs. 9 and 12). 


EQUATIONS 


In the following section all relevant equations (ref. 2) will be given. For the sake of clarity the 
equations are written for a cartesian coordinate system, although the computations are performed with 
the equations transformed into general curvilinear coordinates. 


Mass conservation: 


Momentum conservation: 


Turbulent kinetic energy: 
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Specific dissipation equation: 
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Upper partition stress tensor: 
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In the above notation, u,i is the velocity vector, x i carthesian coordinates, k the turbulent kinetic energy, 
u the specific dissipation rate, p the density, t the time, p the pressure, f %J the sum of the molecular 
and the Reynolds-stress tensor, T x j the upper partition stress tensor and r x j the Reynolds-stress tensor. 
The eddy viscosity vt — Pt/ P is defined as follows: 
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The mean strain rate S x j and the mean rotation tensor are 
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The Reynolds-stress tensor is computed from the upper partition stress tensorT/y: 

rr, 2 c 

r%j = PMj - ~ pebij (9) 

where e is the energy of the eddies in the lower partition. The value of e can be computed by contracting 
equation (9) (note that T u is negative): 

e = k + -Ta ( 10 ) 

The tensor E^j defines the transfer from the upper to the lower partition energies: 

(2 \ 2 / \ 2 
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The tensors P{j and Dfj are production tensors defined by 


The model constants are 
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Ci = 1 +4(1 - e/fc)( 3 / 2 ) (13) 

In case the computations are performed with the two-equation k - u model, the value of £ is £ = 0 
(ref. 15). The Reynolds-stress tensor is then defined as 


T ij — 2 fit S{ j ^ hti'ij (14) 

The author (ref. 11) has recently proposed a modification to the two-equation k-uj model that accounts 
more realistically for the transport of the principal shear-stress component than equation (6). In the 
modified model the eddy-viscosity is defined by 
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where Q is the absolute value of the vorticity and a\ = 0.3. The rationale behind this modification is 
that the equation for the principal shear-stress (-u'v'^j in a boundary layer ( Q = ) can be rewritten 

as follows: 
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In case of an adverse pressure gradient flow the ratio of production to dissipation can be significantly 
larger than one in the outer part of the boundary layer, making the shear-stress larger than a\k, contrary 
to experimental evidence. Equation (15) ensures that < a\k. This approach is very similar 

to the one taken in the Johnson-King model (ref. 13). There a baseline model is supplemented by a 
transport equation for the principal shear-stress that is based on the assumption that —u'v' = a\k in 
the outer part of the boundary layer. In case of the k — u model this assumption can be satisfied by 
applying equation (15). 


For the two-equation models the specific dissipation uj is computed from 
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If the modified version of the k — to model is used, the diffusion constants were changed to 

n* =0.85, (j = 0.6 


The two constants have been optimized to obtain good results for a zero pressure gradient flat plate 
boundary layer. 


The mean flow equations are solved by the INS3D code of S. E. Rogers and D. Kwak (ref. 1) 
which is based on a pseudo-compressibility method. The convective terms are upwind differenced with 
a third-order flux-difference approach. The viscous fluxes are differenced with second-order accurate 
central differences. The linear equations resulting from the first-order backward Euler time differencing 
are solved with an implicit line relaxation scheme. The same methodology used in INS3D has been 
applied to the solution of the turbulence equations. The turbulence equations are solved decoupled 
from the mean flow equations. This means that the turbulence quantities are kept constant during one 
integration step of the mean flow and vice versa. Furthermore, the turbulence equations are solved 
decoupled from one another as scalar equations. This approach is justified by the fact that they are 
coupled mainly through the production and dissipation terms, unlike the much stronger coupling of 
the mean flow equations through the convection and the pressure terms. The advantage of the scalar 
treatment is that it is much easier and faster to invert eight tridiagonal matrices than one 8x8 block 
tridiagonal matrix. 

The main difference in the numerical procedure for the mean flow and the turbulence equations 
is that a total variation diminishing (TVD) scheme is applied to the latter. TVD schemes have been 
developed to prevent oscillations of the mean flow quantities near shock waves in compressible flows. 
Since there are no shock waves in the present incompressible applications, it is not necessary to apply a 
TVD scheme to the mean flow equations. However, the turbulence equations can have weak solutions 
(discontinuity in the derivative) (ref. 14) near the boundary layer edge. A higher order upwind scheme 
can lead to small oscillations in this region. Note that the turbulence quantities approach small values 
at the boundary layer edge and that the oscillations can bring them close to zero. This in turn can cause 
spikes in the ratios of turbulence quantities that appear in the production and dissipation terms. The 
TVD scheme has not only proven helpful in connection with the RSM, but also with the two-equation 
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models, where small undershoots in u can lead to overshoots in Vf Similar effects have been seen with 
k — e models. 

In the following, <j> stands for any of the turbulence quantities of equations (3), (4), and (5) 
and P and D are the corresponding production and destruction terms (divided by J) in that equation. 
Transformation of the equations into a curvilinear coordinate system £ j , i = 1,3 and application of the 
summation convention gives 
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or in short notation: 
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where J is the Jacobian of the coordinate transformation. The effective viscosity z> is equal to v = 
v + (T*"vt for the /c-equation, v — v + for the u equation and z> = 0 for the equations of the upper 
partition stress tensor. 


Spatial Discretization 

The difference scheme for the convection terms is given for the flux / in the £ direction: 
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The flux limited values of the TVD scheme df are defined by the minmod function given in (ref. 6): 
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A value of x — 3 was used for the parameter of the TVD scheme. The minmod function is defined as 
follows: 

minmod[x , y] = sign{x)max [0, min[\x\,ysign(x)]} (26) 

Note that the limited fluxes in equation (25d) are identical to their unlimited values over most of the 
computational domain, and that the limiter is activated only if the fluxes at neighboring points have 
different signs, or are different by a factor of x> The unlimited scheme is third-order accurate for the 
convection terms. 


The viscous terms /* in equation (19) are discretized with standard second-ordercentral differences 
and the production and dissipation terms are simply transformed by applying the chain rule. 


The implicit time discretization of the equations is the same as used in INS3D (ref. 1) for the mean 
flow equations: 
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where n denotes the time level and R the residual of the discretized equations. ^ is an approximate 
Jacobian of the residual R, which accounts only for the first-order terms in the upwind differencing 
ec f (21)) an d the orthogonal terms of the viscous expressions. Forming the exact Jacobian 

is too costly and destroys the tridiagonal structure of the matrices. One of the crucial points in the 
discretization of higher order turbulence models is the linearization of the source terms. Relatively 
good convergence properties have been achieved with the following approximate linearizations: 
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ll was found to be advantageous to use the largest possible value C \ max — 5 for C \ in equation (13) in 
lhe implicit part to enhance diagonal dominance. Note that all the approximations to the implicit part 
of the equations have no effect on the steady-state solution. 


The linear equation system resulting from equation (27) has the following form (ref. 1): 

= R (29) 

where A / is, after appropriate numbering of the equations, a banded matrix of the general form: 


Equation (29) is solved by a Gauss-Seidel line relaxation method: 

A( f’ m+X = & n - m jA i_ l A<p^ kl _ i — m j,k - ] i j 

~ m j,k + 1 1 ,1 ~ m j,k,l+\ A< t > Tk,l+\ (31) 


where m is the iteration level. The tridiagonal system of equation (31) is solved with a LU decomposition 
method. In order to enhance the stability of the solution procedure, the eddy-viscosity formulation is 
maintained in the implicit part of the mean flow equations. 


All boundary conditions are treated implicitly to preserve the stability properties of the implicit 
scheme. At a no-slip condition, the mean flow velocities are set to zero as well as the pressure gradient 
normal to the wall. All turbulent quantities, except u are also set to zero. As pointed out by Wilcox 
(ref. 15), u) satisfies the following equation near the wall: 

W -r-^r as y -> o (32) 

(3y z 

Wilcox recommends to specify this analytical solution for the first few gridpoints away from the wall 
explicitly. The present author found it much easier and as accurate to implement the following boundary 
condition: 

w = io ^(£f at y = 0 (33) 

where Ay is the distance to the next point away from the wall. Equation (33) simulates the boundary 
condition (eq. (32)) without the need of changing the solution at interior points. It should be noted that 
models based on the cj-equation give accurate results if the near wall values of uj are sufficiently large. 
Both equations (32) and (33) satisfy this demand. The results are not sensitive to the factor (10) used 
in equation (33). 
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At a slip surface, the normal pressure gradient is again zero and so is the velocity normal to the 
surface. Furthermore, the gradients normal to the surface of the two velocity components lying in the 
plane of the surface are also set to zero. 

At an inflow and outflow boundary, the mean flow quantities are determined by the method of 
characteristics, This means that for an inflow condition the velocity components are specified and the 
pressure is determined from within the domain, whereas for an outflow condition the pressure is given 
and the velocities are determined by the characteristic relations. The turbulence quantities are specified 
at inflow boundaries; at outflow boundaries, a zero gradient is assumed. 

Two of the computed flowfields have a rotational symmetry. In these cases, the gradients of all 
quantities in the circumferential direction are set to zero. 

RESULTS 


The first set of computations has been performed for a flat plate zeropressure gradient boundary 
layer. The number of gridpoints was 60 x 3 x 80. Although the flow is two-dimensional, it is necessary 
to introduce three planes in the lateral direction since the code solves only three-dimensional flow 
problems. Resolution studies with a 40 x 3 x 80 and a 90 x 3 x 90 grid show that the solutions are grid 
independent. All grid refinement studies have been done with the k — u model for economic reasons, 
but are believed to be valid also for the multiscale model because of the strong similarities between 
the equations. Figure 1 shows the results of the computations. Obviously, all models give very similar 
results, both in the cj -distribution and in the law-of-the-wall plot. It should be mentioned at this point, 
that the results computed with the k — cj model have a strong dependence on the freestream values of 
cj, as shown in reference 16. In order to ensure unambigious solutions, the proper inflow values for 
cj given in reference 16 have been specified. It should be noted that the freestream dependency of the 
model is not removed by the modifications to the eddy-viscosity given in equation (15). 




Figure 1. Flat plate boundary layer, (a) Skin friction, (b) law-of-the-wall. 
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Adverse Pressure Gradient Flow 

The second test case is an adverse pressure gradient boundary layer, developing in the axial direction 
of a circular cylinder, as reported by D. Driver (refs. 7 and 8). The pressure gradient is strong enough 
to cause the boundary layer to separate. The computations have been carried out with a 60 x 3 x 60 
grid. The grid independence was ensured by a 100 x 3 x 100 gridpoint computation with the k — u 
model. 

This flowfield has been computed before by the present author (ref. 1 1 ) with a number of different 
eddy-viscosity models. The best results in that comparison have been obtained with the Johnson-King 
(JK) model (ref. 13). The major difference between the JK model and equilibrium algebraic eddy- 
viscosity models is the inclusion of the transport effect of the principal turbulent shear-stress component. 
One would expect that the RSM and the modified k — to model should be able to improve the results 
of the baseline k — uj model, since both of them account properly for this transport effect. Note that 
the modified k — uj model is based on the same assumption as the JK model, that the shear-stress is 
proportional to the kinetic energy in the outer part of the boundary layer. 

Figure 2 shows the computed c p distributions in comparison with the experimental data. It is 
obvious that both the MS model and the modified k — cj model give superior results over the baseline 
k — uj model. The predictions of the two former models are virtually identical and in excellent agreement 
with the experimental data. Keeping in mind that the c p distribution is the single most important quantity 
in aerodynamic computations, the distribution shown in figure 2 gives some hope that properly designed 
two-equation models and full RSM can give accurate predictions for these flows. A comparison of 
computed and measured cy distributions is shown in figure 3. There are only small differences in 
the predictions of the k — uj and the MS models, whereas the modified model predicts the minimum 
wall-shear-stress more accurately, but overpredicts the downstream extent of the separation zone. 

The differences seen in figures 2 and 3 are consistent with the velocity profiles depicted in figure 4. 
The MS and the modified k — uj model predict higher displacement thicknesses than the k — uj model 



x/D 

Figure 2. c p distribution for adverse pressure gradient flow CS0. 
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x/D 

Figure 3. Cf distribution for adverse pressure gradient flow CSO. 




Figure 4. Velocity profiles for adverse pressure gradient flow CSO. 


and therefore a superior pressure distribution. Near the wall the modified model is in better agreement 
with the experimental shape of the profiles and therefore gives lower cj values. 


10 







Figure 5. Shear-stress profiles for adverse pressure gradient flow CSO. 

The main reasons for the differences in the predictions can be seen in figure 5, showing two 
turbulent shear-stress profiles at different axial locations. It was found in all computations of this 
flowfield that the c p distribution can be predicted accurately if the shear stresses up to separation are 
computed correctly. Note that the stresses predicted by the MS model and the modified k — c o model 
are in close agreement with experiments at x/D = 0.363 and that the differences in the predictions at 
x/D = 1 .633 are not reflected in the c p distribution (fig. 3). The reason for the overprediction of the c p 
values by the k — u model lies in the overprediction of the turbulent shear-stress levels up to separation 
(fig. 5) and is an immediate consequence of equation (6). 

Figure 6 shows the distribution of the normal stresses, as predicted by the MS model. Note that in 
the present case, v ' 2 is the stress component in circumferential and w /2 the one in wall normal direction. 
The normal stresses are obviously predicted quite well, although slightly too low at x/D = 0.363. 

The above results confirm the theory motivating the JK model that the main additional complexity 
in the prediction of adverse pressure gradient flows comes from the importance of the transport terms in 
the equation for the principal shear-stress. The comparison of the results from the MS and the modified 
k — uj model shows that a properly designed two equation model can give as accurate results for this 
kind of flow as a full Reynolds-stress model. 
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Figure 6 . Normal-stress profiles for adverse pressure gradient flow CSO. 


Spinning Cylinder Flow 

The next problem that was considered is the flow along a circular cylinder with a spinning upstream 
{x/ D < 0) and a stationary downstream (x/D > 0) section, as reported by D. Driver and J. P. Johnston 
(refs. 8 and 9). The circumferential surface velocity of the spinning part is equal to the incoming 
freestream velocity Uqq in axial direction. 

The flow is collateral in the upstream (rotating) part of the cylinder, as noted by Driver et al., which 
means that the velocity profiles are two-dimensional for an observer moving with the wall. Therefore 
no three-dimensional effects have to be expected in this region. However, due to strong curvature 
effects the computation of the upstream flow is by no means trivial. As the flow relaxes back to a 
purely two-dimensional boundary layer along the stationary part of the cylinder, the three-dimensional, 
curvature and transport effects are of equal importance, which creates a very challenging problem for 
all turbulence models. 

The computations have been performed on a 90 x 3 x 90 grid. Solutions on a 60 x 3 x 60 grid are 
only slightly different, confirming the grid independence of the present results. 

Figure 7 shows two profiles of the axial velocity at different axial stations. The profile on the 
left is located on the spinning part of the cylinder and the profile on the right in the three-dimensional 
interaction region on the stationary part. Figure 8 depicts the circumferential velocity component at 
the same two stations. Obviously, the profiles in the spinning part of the cylinder x < 0 are not 
computed accurately by any of the models. The near wall gradients of both the u and the v component 
at -0 = —0.091 are severely underpredicted. The reason for this discrepancy is that the models do not 
pick up the destabilizing effect of the curvature and therefore underpredict the turbulent shear-stresses, 
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Figure 7. Profiles of axial velocity for spinning cylinder flow AS1. 
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Figure 8. Profiles of circumferential velocity for spinning cylinder flow AS1. 
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as seen in figure 9. The MS model gives somewhat higher shear-stress levels than the two-equation 
models, but failed to meet expectations (refs. 17 and 18) of predicting curvature effects properly. 

With regard to the downstream station & = 0.726, the results of the MS model are also not 
significantly better than the results of the two-equation models. The shape of the circumferential 
velocity profile is computed in closer agreement with the experiment (fig. 8), but the absolute values 
are also significantly higher near the wall. The more realistic shape of the v profile is a result of the 
improved prediction of the v'w' profiles (fig. 9). 

Figure 10 shows a comparison of computed and measured turbulent kinetic energy profiles. Note 
that the results of the MS model and the two-equation models are very close to each other and signif- 
icantly smaller than the experimental values, again due to the underprediction of the curvature effects. 
Figure 1 1 depicts the results of the MS model in more detail by showing the individual normal stresses 
compared to the experiments. 

Figure 12 shows the Cf -distribution along the stationary part of the cylinder. The wall-shear-stress 
component in axial direction cj x is underpredicted and the component in circumferential direction Cf y 
is overpredicted, which is in agreement with the velocity profiles (fig. 7). 

Figure 13 shows a map of the second (II) and third (III) invariant of the anisotropy tensor b{j at 
two downstream stations. The definitions of the anisotropy tensor and its invariants are 




Figure 9. Profiles of axial iu'w’) and circumferential {v'w') shear-stress profiles for spinning cylinder 
flow AS1. 
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Figure 10. Profiles of turbulent kinetic energy for spinning cylinder flow AS1. 
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Figure 11. Profiles of normal turbulent stresses, computed from the multiscale model for spinning 
cylinder flow AS1. 
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Figure 12. Axial ( Cf x ) and circumferential (cf y ) wall-shear-stress distributions for spinning cylinder 
flow AS1. 



Figure 13. Anisotropy invariant map, computed from the multiscale model for spinning cylinder flow 
AS 1 . 
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( 34 ) 


^ u iUj ^ 

bij ~ 2k 

U = -bijbji/2 , /// = bijb jk b ki /3 (35) 

The map was introduced by Lumley (ref. 1 9) to characterize the anisotropy of turbulent flows in a coor- 
dinate invariant fashion. Each of the boundaries of the map corresponds to a different asymptotic state 
of turbulence: the origin is the isotropic state, the upper boundary corresponds to all two-dimensional 
states, the left boundaiy represents an axisymmetric contraction and the right boundary an axisymmetric 
expansion. It is obvious from figure 13 that the predictions of the MS model cover a significantly smaller 
part of the map than the experimental data. Although the flowfield is highly three-dimensional and has 
strong curvature effects, the predictions stay very close to their design point for a zero-pressure-gradient 
flat plate boundaiy layer. The inability of the model to predict the anisotropic behavior of the flow is 
suspected to be a shortcoming of the Launder-Reece-Rodi pressure-strain model (ref. 3) underlying the 
MS model. Results of boundary layer computations, based on the Launder-Reece and Rodi model as 
reported by Driver (ref. 8), are in very good agreement with the present results, giving further support 
to this opinion. 


The flowfield discussed in this section has been investigated experimentally by S. D. Anderson and 
J. K. Eaton (ref. 10). The experimental setup is shown in figure 14. All dimensions are in centimeters. 
A two dimensional boundary layer enters the test section and is deflected by a 90° wedge. Due to 
the strong lateral pressure gradient created by the wedge, highly skewed three-dimensional boundary 
layer profiles develop downstream. Measurements of the mean velocity and of turbulent correlations 
are available along the centerline, and along a selected streamline, as shown in figure 14 (experimental 
ports). In addition to the profile information, lateral static pressure distributions have been measured 
at seven x-locations outside the boundary layer. In the following, the origin of the coordinate system 
is located at the tip of the wedge, with the x-coordinate pointing in the direction of the symmetry line 
(positive x in direction of the incoming flow), the y-coordinate is in the lateral and the z-coordinate in 
surface normal direction. 

The geometry of the test section is symmetric with respect to both the y- and the z-direction, so 
that in principle, only one quarter of the flowfield needs to be computed. Unfortunately, a small amount 
of suction was applied in the experiment to the incoming boundary layer on the upper wall, so that the 
upper and lower incoming boundary layers were slightly different. 

Three different sets of computations have been performed with the k — lj model in order to test 
the grid dependence and the influence of the asymmetric inflow conditions. First a grid refinement 
study was performed computing only one quarter of the flowfield. It was found that the computations 
performed on a 35 x 35 x 35 grid and on a 50 x 50 x 50 grid are in very close agreement with each 
other. Figure 15 shows a comparison of streamwise and the crossflow mean velocity profiles at the last 
downstream station of the experiments. Differences further upstream are even smaller. An additional 
computation was performed on a 40 x 40 x 80 grid, including the upper boundary layer. Results of that 
computation are also included in figure 15. The inclusion of the upper half of the flowfield obviously 
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Figure 14. Experimental setup for 90° wedge flow (from ref. 10). 








does not lead to significant changes in the solution. The following results are therefore all for the 
50 x 50 x 50 grid with a symmetry condition in the middle of the tunnel. 

Figure 16 shows a comparison of computed and experimental pressure distributions for the three 
different models. Although the differences are not dramatic, the MS- and the modified k — u model are 
in better agreement with the experiments than the original k — uj model, especially in the downstream 
region. 

Figure 17 shows the velocity profiles along the centerline of the channel. Differences between the 
models appear only at the station closest to the tip of the wedge. The k — u model does not react 
sufficiently to the adverse pressure gradient and predicts too full a profile shape. The modified k — uj 
model is separated at this station, and the MS model is close to separation. Although the MS model is 
apparently in better agreement with the data, no final conclusion can be drawn from figure 17, because 
the three hole probe used in the experiments cannot measure separated profiles. Unfortunatelly, no 
oil-flow pictures are available for this experiment. 

A comparison of turbulent shear stress profiles is given in figure 18. Experimental data are only 
available for the first three stations. The larger symbols at z = 0 represent the measured wall-shear 
stress. The results are very similar to those for the adverse pressure gradient flow csO (fig. 5). The 
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Figure 16. C p distribution for 3D wedge flow. 
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Figure 18. Shear-stress profiles along centerline for 3D wedge flow. 
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k — lu model overpredicts the shear stress levels, whereas the other two models are in close agreement 
with the experiments which in turn leads to more retarded velocity profiles. 

The present computation focuses on the performance of the models in the region of three-dimensional 
boundary layer profiles away from the centerline. Figure 19 shows four streamwise velocity profiles 
along the selected streamline. The local coordinate system for these plots is aligned with the experi- 
mental freestream direction. Figure 20 shows the corresponding crossflow profiles. The modified k — to 
model is obviously closest to the experiments. It predicts the retardation of the streamwise velocity 
component near the wall, and gives the correct amount of crossflow at the last downstream station. The 
discrepancies between the computations and the experimental profiles in the upstream portion of the 
flow might be explained by the failure of all of the models to predict the reduction of the turbulent 
shear stress shown in figure 21 due to the three-dimensionality of the mean-flow profiles. However, 
both the MS and the modified k — to model predict significantly lower shear-stress values than the 
original k — u model and are therefore closer to the experimental data. Note that experimental stress 
profiles are only available for the first three stations. The close agreement between the former two 
models in the prediction of the turbulent shear stresses shows that the advantage of the MS model 
over the baseline two-equation model is again its ability to account for transport effects. This is also 
demonstrated by figure 22, showing the shear stresses in the crossflow direction. Again, the MS and the 
modified two-equation model are in very close agreement with each other and somewhat closer to the 
experiment. Obviously all models fail to predict the crossflow component of the shear-stress correctly 
and are therefore not able to predict the correct direction of the shear-stress vector. This has to be 
expected from the eddy-viscosity models which assume that the stresses are proportional to the strain 
rate. The MS model in turn solves a separate equation for every stress component and has therefore the 




Figure 19. Streamwise velocity profiles along selected streamline for 3D wedge flow. 
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Figure 21. Streamwise shear-stress profiles along selected streamline for 3D wedge flow. 
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Figure 22. Crossflow shear-stress profiles along selected streamline for 3D wedge flow 


potential to give more accurate results, but in practice it is not more successful than the two-equation 
models. 

The above finding that the inclusion of the transport effects improves the prediction of pressure- 
driven three-dimensional boundary layer computations are in agreement with Kavasaoglu et al. (ref. 20). 
In their paper the pressure-driven boundary layer flow of Van den Berg and Elsenaar (ref. 21) was 
computed with different turbulence models. The Johnson-King model (ref. 13) improved the predictions 
considerably when compared with models that did not account for the transport of the principal shear 
stress. 


A numerical method has been presented to solve a full Reynolds-stress model in combination with 
the three-dimensional Navier-Stokes equations. The equations for the Reynolds-stresses are solved 
decoupled from the mean flow equations with a higher order upwind TVD scheme. The numerical 
stability of the algorithm is enhanced by retaining the eddy-viscosity formulation and treating the rest 
of the stresses as explicit source terms. 

Computations have been performed for four different flowfields, and the numerical results were 
compared with solutions based on the k~ uj model in its original form and in a modified version recently 
proposed by the author. The modified version accounts for the transport effect of the principal shear 
stress by limiting the eddy-viscosity in a way that keeps the principal shear stress smaller or equal to 
03k. 
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As expected, all three models give accurate predictions for the zero-pressure-gradient flat plate 
boundary layer. 


The computations for the separated adverse pressure gradient flow show a decisive advantage of 
the MS and the modified k — u model due to their ability to properly account for the transport of the 
shear stresses. 

For the spinning cylinder flow, all three models produce rather unsatisfying results. None of the 
models is able to predict the destabilizing effect of the curvature on the turbulence. While this is to be 
expected from the eddy- viscosity models, it is also true for the Reynolds-stress model. Obviously, the 
correct treatment of the production terms is not sufficient to predict this complicated phenomenon. 

The last set of computations was performed for a three-dimensional wedge flow. In this flow a 
strong lateral pressure gradient produces highly skewed velocity profiles. The main challenge of this 
type of flow is that the shear-stress direction is no longer aligned with the strain-rate direction, and that 
the skewing of the velocity profiles reduces the shear-stress levels considerably. Both effects are beyond 
the scope of the eddy-viscosity models, but are also not predicted by the MS model. The success of the 
modified k — uj model in predicting the correct amount of crossflow shows, however, that the reduction 
of the shear stress levels due to the inclusion of the transport terms is at least as important an effect, 
and gives some hope that this type of flow can be predicted with sufficient accuracy by a comparatively 
simple model. 

In conclusion, the full Reynolds-stress model did not show an advantage over an improved two- 
equation model for the flows considered in this paper. Its main advantage over standard two-equation 
models seems to be the ability to account for the transport terms of the turbulent stresses. Obviously, 
this can also be achieved within the framework of two-equation models. The additional complexity of 
solving a full Reynolds-stress model does not seem to be justified for the types of flows considered 
here, unless better models can be developed. 
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